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Improved Multilevel Monte Carlo Methods for Finite Volume 
Discretisations of Darcy Flow in Randomly Layered Media 

M. Park* and A.L. Teckentrup^ 


Abstract 

We consider the application of multilevel Monte Carlo methods to steady state Darcy flow in a 
random porous medium, described mathematically by elliptic partial differential equations with random 
coefficients. The levels in the multilevel estimator are deflned by finite volume discretisations of the 
governing equations with different mesh parameters. To simulate different layers in the subsurface, the 
permeability is modelled as a piecewise constant or piecewise spatially correlated random field, including 
the possibility of piecewise log-normal random fields. The location of the layers is assumed unknown, 
and modelled by a random process. We prove new convergence results of the spatial discretisation 
error required to quantify the mean square error of the multilevel estimator, and provide an optimal 
implementation of the method based on algebraic multigrid methods and a novel variance reduction 
technique termed Coarse Grid Variates. 


1 Introduction 

Mathematical models of physical processes are frequently used for simulation. The parameters appearing 
in these models are often subject to uncertainty, due to for example lack of measurements or data, and the 
result of the simulation hence also becomes uncertain. In this work, we are interested in the simulation of a 
steady state groundwater flow, governed by Darcy’s law, in a porous medium of which the permeability is 
not fully known. Assigning a suitable probability distribution to the permeability, the goal of the simulations 
will be to compute moments of quantities of interest related to the resulting pressure head and Darcy flux. 

In recent years, multilevel Monte Carlo (MLMC) methods have frequently been applied and analysed in 
the context of partial differential equations with random coefficients 0 [231 [Ml E na Ian 1 [201 [M] . Originally 
introduced by Heinrich [22] in the context of parameter dependent integral equations and by Giles m in the 
context of stochastic differential equations arising in mathematical finance, multilevel Monte Carlo methods 
present a significant computational saving compared to standard Monte Carlo methods. Exploiting the 
linearity of expectation, multilevel Monte Carlo methods combine a large number of samples of the quantity 
of interest at low spatial resolution with a small number of samples at finer spatial resolutions, to provide 
an accurate estimate of the quantity of interest at fine spatial resolution at low computational cost. 

The purpose of this paper is to extend the theoretical convergence results available for MLMC methods, 
as well as improve on their practical performance through a novel variance reduction technique, with a 
particular focus on aspects relevant to the groundwater flow application. On the theoretical side, we firstly 
provide new convergence results on the spatial discretisation error, required to bound the mean square error 
of the MLMC estimator, in the case of finite volume discretisations. In applications such as groundwater 
flow modelling, finite volume methods are often preferred over methods such as standard finite elements due 
to local mass conservation . Secondly, we extend the range of layered permeability models considered in 
[M] by allowing also the location of the layers to be uncertain. For the sake of generality, throughout the 
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analysis we do not assume uniform coercivity or boundedness of the permeability in terms of the random 
parameter, but instead follow the more general framework suggested in [2123]. 

On the implementation side, we propose a novel variance reduction technique, termed Coarse Grid 
Variates, to further lower the computational cost of MLMC estimators. The technique is designed to work 
for stationary models of the permeability. Assuming that for each sample of the permeability, we compute its 
values at the cell centres of the finite volume grid, we argue that samples of the permeability on a coarse grid 
can be extracted from the information contained in a sample of the permeability on a finer grid. The exact 
number of coarse grid samples that can be extracted from a single fine grid sample depends on the ratio of 
degrees of freedom between the two grids. Using this observation, together with an averaging procedure, we 
formulate a new MLMC estimator which is provably unbiased and in numerical simulations shows a variance 
reduction of up to 2 orders of magnitude at less than twice the computational cost of the standard MLMC 
estimator. 

The outline of the remainder of this paper is as follows. In section|2 we present the mathematical model 
of interest, together with any assumptions that we make on its components. Section [3] starts with a brief 
overview of finite volume methods and MLMC estimators, before we in section 13.31 provide a convergence 
analysis of MLMC estimators based on finite volume discretisations. The theoretical results are illustrated 
in two and three spatial dimensions in section 13.41 Section |2 is devoted to the derivation and numerical 
simulation of the Coarse Grid Variates technique. Finally, section |2 provides some conclusions. 


2 Problem setting 

The study of groundwater flow is well established, and there is general scientific consensus that in many 
situations Darcy’s law can be expected to lead to an accurate description of the flow [12 min- The 
classical equations governing (steady state) single phase subsurface flow consist of Darcy’s law coupled with 
an incompressibility condition: 

q + KVp = g and div q = 0, in D d=l,2,3, (2.1) 

subject to suitable boundary conditions. In physical terms, p denotes the pressure head of the fluid, K is 
the permeability tensor, q is the filtration velocity (or Darcy flux) and g are the source terms. 

A typical approach to incorporating the uncertainty in p and q is to model the permeability tensor as 
a random field K = K(x, w) on D x D, with (D,-F, P) a probability space [TT] [T^. The model (12.11) then 
becomes a system of partial differential equations (PDFs) with random coefficients, which can be written in 
second order form as 

— div(K(x, a;)Vp(x, w)) =/(x), in D, (2.2) 

with / = —div g. The solution p itself is also a random field on D x D. For simplicity, we assume 
that the boundary conditions and the sources g are deterministic, and restrict ourselves to convex polygo¬ 
nal/polyhedral domains D. For the more general case, we refer the reader to [231122 . 

In this general form, solving (ESI) is extremely challenging computationally, and in practice it is therefore 
common to use relatively simple models for K that capture the most important features of the subsurface 
geometry and are as faithful as possible to the measured data. For simplicity, we restrict our attention to 
scalar valued models of the permeability, i.e. we replace the tensor K by a scalar valued function k. Tensor 
valued coefficients are considered in [24) . In this paper, we are particularly interested in discontinuous models 
of the permeability, which are of particular practical interest due to their ability to model the different layers 
present in the subsurface. To this end, we assume that the computational domain D is partitioned into m 
disjoint convex polygonal subdomains with the permeability k continuous on each subdomain Di, 

i = 1 ,..., m. 

In applications the exact location of the different layers of the subsurface is often not known exactly, and 
hence also uncertain. We incorporate this additional uncertainty into our model by allowing the partitioning 
{DilEi to ^Iso be random, independent of k, leading to a random partitioning {Di(w)}Ei) for a; G D. As 
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before, we assume that for each realisation of the random partitioning, i.e. for a fixed w, the m subdomains 
are disjoint and convex polygonal. 

It remains to specify the model of the permeability employed in the subdomains. For each subdomain, 
denote the permeability k restricted to Di hy ki. A model that has been studied extensively in subsurface 
flow applications is the log-normal distribution, which allows the permeability to vary over many orders of 
magnitude and guarantees that the permeability takes positive values for almost all realisations. We will in 
particular work with the following two models: 

Example 1. (Piecewise constant model) In each subdomain Di, we model the permeability as a log-normal 
random variable: ki{x.,uj) = exp{Zi{uj)), where Zi ~ N{p,i,af), i = 1,... ,m. The mean and variance of the 
random variable are allowed to vary between the subdomains. As a function of the spatial variable x, each 
ki is constant. 

Example 2. (Piecewise correlated field model) In each subdomain Di, we model the permeability as a log¬ 
normal spatially correlated random field: fci(x, w) = exp(gi(x, w)), where gi is a stationary Gaussian random 
field with constant mean /ii(x) = fii and an exponential two point covariance function: 

C{x,y) = E[(g,(x) - ^J.i){g^{y) - tn)] = exp(-||x - y||r/A*). (2.3) 

Here, || • ||r denotes the usual G norm on and typically r = 1,2. The parameter Xi is known as the 
correlation length and af as the variance of the Gaussian field gi. It can he shown that in this case ki, as 
a function of the spatial variable x, is Holder continuous with exponent less than 1/2, ki{-,uj) G C*{Di), for 
any t < 1/2. 

Other models of the permeability are of course possible, and the results of this paper are readily applicable 
in a variety of situations. In particular, the analysis in this paper holds also in the case where the covariance 
function C{x, y) in Example [2] is replaced by another member of the Matern class of covariances with 
smoothness parameter > 1/2 [4]. 


3 Finite volume methods and multilevel Monte Carlo sampling 

In this section, we will describe and analyse the numerical methods employed in this paper. We will start 
with a description of the finite volume method used for the spatial discretisation in section 13.11 before we 
briefly recall the multilevel Monte Carlo method in section lX^ and finally in sectionprove the convergence 
of the multilevel Monte Carlo method for the finite volume method considered in section EH 


3.1 Finite volume discretisation 

The starting point of finite volume discretisations is the second order formulation (EU. One then chooses a 
non-overlapping partitioning of the domain D into boxes (or volumes) Bh, where h denotes the mesh width 
of the partition. Integrating equation (EB over each box B G Bh leads to a set of algebraic equations 


— f div(fc(x, a;)Vp(x, w)) dx = j /(x)dx, \/B 

J B J B 


GBh. 


(3.1) 


The volume integral on the left hand side is transformed into a boundary integral using the Divergence 
Theorem: 

(3.2) 


/ /c(x, w)Vp(x, w) • nds = / /(x) dx, 

JOB J B 


where n denotes the unit outward normal and dB denotes the boundary of the box B. The specific finite 
volume scheme is now determined by the choice of volumes Bh, as well as how the integrals in (13.21) are 
computed (exactly or by quadrature). 

We will in this paper consider cell-centred finite volume methods on uniform rectangular meshes. For 
illustrative purposes, let us describe this discretisation in more detail in the particular case where the 
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computational domain is the two-dimensional unit square, D — (0,1)^, and the boxes Bh are squares. The 
cases of rectangular domains, and one or three spatial dimensions, are treated analogously. 

We start by subdividing [0,1]^ uniformly into a mesh of m x m square cells and denote by Bi^j the cell 
(^’ m) ^ m)’ = 1,... ,TO, and by its centre. We have Bh = with m = 1/h. 

Let kij and fij be the values of k and / at Xij, respectively, and denote by pJJ our approximation to p at 
Xij . We will approximate the right hand side of (EH) by the midpoint rule. 



/(x)dx 




To approximate the left hand side of E3D we consider separately each edge of dBij. The contribution from 
the edge between Bij and i3i+ij is again approximated by the midpoint rule. A simple approximation to 
k on the edge is its value at the midpoint, k^j^i y As approximation to the gradient Vp • n on the edge we 
use the central finite difference (pf^i j The contributions from the other edges are approximated 

similarly, leading to the following form of the (i,j)-th equation: 


- - K+i,JP^+l,3 - p^ J+l = ko! 


(3.3) 


where j + j ^ Neumann boundary condition, i.e. a prescribed flux 

—kVp ■ n = i/jat, on any part of the outer boundary of (0,1)^ is straightforward to incorporate. The 
respective flux term on the left hand side of (13.3p is simple replaced by V'at evaluated at the midpoint of the 
edge. To enforce a Dirichlet boundary condition, i.e. a prescribed pressure p = 1l}D^ we simple replace the 
central difference by a one-sided difference. 

An alternative approximation to (13.3p . which is often used in subsurface flow applications, is derived 
using the harmonic average of kij and fci+ij as the value of k on the edge between Bij and Bi+ij. 

The analysis in section 13.31 can be applied also in this case, leading to the same convergence rates as in 
Theorem 13.21 

The linear system of equations arising from the approximation (13.31) takes the standard five-point stencil 
form. It is sparse and highly ill-conditioned, but it can be solved efficiently and robustly with algebraic 
multigrid (AMG) methods. In fact, we will see in section 13.51 that an iterative solver based on an AMG 
preconditioned Gonjugate Gradient (GG) method scales optimally even in 3 spatial dimensions. 


3.2 Multilevel Monte Carlo methods 

We now briefly review the ideas of the multilevel Monte Garlo (MLMG) technique. For more details, we 
refer the reader to [13 [5]. 

Suppose we are interested in finding the expected value of a functional Q — Q{p), where p is the solution 
to the Darcy flow equation (12.21) . Examples of functionals Q of interest include the value of the pressure head 
p or the Darcy flux —kVp at a particular point in the computational domain D, or the outflow over parts of 
the boundary. Since p can not be computed exactly, we in practice use a finite volume approximation of Q, 
denoted by Qh := 0(Ph^)- 

The main idea behind the MLMC technique is now as follows. Gonsider simulations with different mesh 
widths hg, chosen such that 

hi = for £ = 0,1, ...,L, (3.4) 

where s is a positive integer. In contrast to the standard Monte Carlo (MC) approach, which only uses 
samples of Qh generated on the finest level L, samples on all grid levels £ = 0,..., L are taken into account 
in MLMC to estimate statistical moments of solution. Using the linearity of expectation, and denoting 
Qt := Qhi, we have 

L 

E[Ql] =E[Qo]+^E[g^-Qf_i]. (3.5) 

1=1 
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Each of the expectations on the right hand side of (13.51) is now estimated independently using a Monte Carlo 
estimator, resulting in the MLMC estimator 


L Ni 

t=0 ^ ^ i=l 

where for simplicity we have set Q-i = 0. The number of Monte Carlo samples Ni, on each level is chosen 
such that the overall variance of the multilevel estimator is minimized for a fixed computational cost. It is 
important to note that the quantity Qty — Qi-i,i in (13.61) is computed from two discrete approximations 
with different mesh widths, but the same random sample 

In order to quantify the accuracy of the multilevel estimator (13.61) . we consider the mean square error 
MSE(Q, Q) of the estimator Q as an estimator of Q- 

MSE(Q, Q) = V[g] + (E[Q] - Qf. 


Using the unbiasedness of Monte Carlo simulations, together with the fact that the L + 1 individual Monte 
Carlo estimators in (13.61) are independent, the mean square error of the MLMC estimator is 


MSE(Qr.Q) + (3.7) 

£=0 ^ 

To achieve a mean square error MSE(Q^'^,(5) at the tolerance level e^, we evenly distribute between 
the two terms on the right hand side of (1X71) . We furthermore denote the cost to compute one sample 
Qi,i — by Cl, which includes both the cost of producing the sample of k and solving the corresponding 

finite volume equations. We have the following results on the computational cost of the MC and MLMC 
estimators to achieve a mean square error of e^. 

Theorem 3.1. Suppose there exist positive constants a, > 0 such that a > bmin(/3,7) 


{Ml) \E[Qi-Q]\<Cc,h'^, 
{M2) Y[Qi-Qe-i]<Cph^,, 
{MS) Cl < C^hTg^^. 


Then for any e < e there exists a positive constant a value 

MSE(Q“l,Q) < e^, and 


( ifP 

C,{Qf^) = { C“Lg-2(loge)2, *//3 

^ (7MLg-2-(7-/3)/a^ ^ 


L and a sequence {fV^}|Uo 


< 7, 
= 7. 
> 7. 


such that 


whereas 


C^{Q^^) = C'“Cg-2-7/« 


for some positive constant . 


A proof of the above Theorem can be found in [ 6 ]. Further reductions in the computational cost of the 
MLMC estimator may be possible by using an optimal, uneven splitting of the total error between the two 
error contributions in (1X71) [ 8 ], and by using optimised mesh hierarchies {hi}^^^ [21]. The rates a, j3 and 
7 are application dependent. The rates a and /3 generally depend on the spatial regularity properties of Q 
and the numerical method used for the approximation Qg. The rate 7 , on the other hand, depends on the 
method of sampling and the method used to solve the linear system of equations for each sample. In the 
best case, we have 7 approximately equal to d, the spatial dimension of the problem. 
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3.3 Finite volume error analysis 

This section is devoted to proving assumptions Ml and M2 of the complexity theorem for finite volume 
discretisations of the model problem (12.11) with the permeability as described in Examples [T] and [5] This 
will be done by showing that the cell-centred finite volume discretisation described in section 13.11 in fact 
is equivalent to a finite element discretisation with a particular quadrature scheme used to assemble the 
stiffness matrix. The convergence of the finite volume discretisation then follows from the convergence of 
finite element discretisations proven in mESlEl]. We would like to point out here that a wider range of 
finite volume methods can in fact be analysed by comparison to a related finite element method, see e.g. 

m- 

For ease of presentation, we again consider the case D = (0,1)^ in detail, and restrict our attention to 
homogeneous Dirichlet conditions tpo = 0 on the entire boundary dD. A standard finite element approxima¬ 
tion of model problem (EB) starts with the weak formulation of (EB, obtained by multiplying the equation 
by a test function v G Hq{D), integrating over the computational domain D and applying Green’s formula: 
find p{-,uj) G Hq{D) such that 

f fc(a:,a;)Vp(x, w) • Vu(x) dx = f /(x)v(x)dx, \/v G Hq{D). (3.8) 

J D J D 

Here, Hq{D) is the usual Sobolev space of functions with square integrable weak derivatives that vanish on 
the boundary: 

Hq{D) := {v : / |up -I- |Vup dx < cx) and u|aD = 0}. 

JD 

The finite element approximation to (j3.8L denoted by , is then defined by p^^{-,oj) G 14 and 


[ /c(x,w)Vp™(x,a;) • Vz;/i(x) dx = [ /(x)uft(x)dx, 'ivhGVu, 

JD JD 


(3.9) 


where 14 is a suitably chosen finite dimensional subspace of Hq{D). For our purposes, we shall choose 
14 to be a space of continuous, piecewise bilinear functions on D that vanish on the boundary dD. To 
facilitate the comparison with the finite volume discretisation described in section 13.11 we choose the mesh 
for the finite element discretisation such that the degrees of freedom in the finite element method coincide 
with the degrees of freedom of the finite volume discretisation, which are the centres x^^ of the boxes 
^ degrees of freedom of the standard, piecewise bilinear finite element 

method are at the vertices of the mesh. This means that for a fixed m = 1/ft,, the finite element mesh is 

given by Bu = where with the nodes ?/° = 0,?/™ = 1 and = [i — l/2)ft, the square elements 

are given by Bij = (j/*, j/*+^) x ( 2 / 42 /-’+^). 

Since we impose Dirichlet boundary conditions on dD, only the interior nodes of the finite element mesh 
are considered as degrees of freedom. Note that these nodes are located exactly at the cell centres x^^-, for 
1 4 4 w. To solve equation (13.9|) . we now choose a basis for the piecewise bilinear finite element space 

14 . We will use the well-known hat functions 4>i,j, whose support is contained in the elements neighbouring 
the node x^^-. In particular, for any point (xi,X 2 ) G [0,1]^, we have 

if (xi,a;2) G Hi-ij-i, 
if {xi,X2) G 
if (xi,a;2) G 
if (xi,X2) G Bi^j, 
elsewhere. 

The finite element solution p™ can be expanded in the basis 4>i,j as p^ = YhiKi j<mPl^ where the 

coefficients pf® will depend the particular realisation, i.e. on w. For ease of presentation, we will from now 


{xi,X2) = < 


(xi-y' ^)(y^-X2) 

(yi-yi-l^l^yj-yj-l) 

(y'-xi)(x2-v^ 

(y' -Xi)(y^ -X2) 

0 
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on drop the dependence of fc, p and p^ on w. By choosing Vh = (j>k,i in (13.91) . for 1 < k,l < m, we obtain a 
linear system of equations for the unknown coefficients pff- 

pff [ k{x)V(l)ij{x)-\/(l)k,i{x)dx= [ f{x)(l)k,i{x)dx, for 1 < i,j,k,l < m. (3.10) 

Jd J d 

Due to the local support of the basis functions, the integrals appearing on the left hand side of (13.101) are zero 
if — fc| > 1 or |j — > 1. We will now devise a quadrature scheme to compute these integrals such that the 

linear equations (13.101) are the same as the linear equations (13.31) defining the finite volume approximation 
denote the corresponding solution of the finite element equations with quadrature by pf^^, 
with coefficients pf^^ ■ 

Let us start with the case k = i + 1,1 = j. Due to the local support of the basis functions, the integral 
over D on the left hand side of (13.101) reduces to an integral over Bij-i U Bij. Expanding the dot product, 
this integral splits into the sum of two integrals. We approximate the first integral, which involves derivatives 
with respect to the first coordinate direction xi, by the midpoint rule in xi and the trapezoidal rule in X 2 ' 

f d<pi+l,j J , / (i+1/2) (j), / (z+1/2) 0 + l)A 

The integral over Bij-i is approximated in the same way, with the integrand evaluated at the two points 
{xi'^^^‘^\x 2 ~^'^) and The second integral, involving derivatives with respect to X 2 , is similarly 

approximated by the midpoint rule in X 2 and the trapezoidal rule in xi. Explicit computations now show 
that this leads to the same set of equations as (13.31) (see, for example [SJ Exercise 4.1.8] and [31 §3.3]). 

The computations above are easily extended to general rectangular domains D and to three (or one) 
spatial dimensions. The quadrature scheme which makes the approximate finite element and the finite 
volume solution equivalent, is the one which uses the midpoint rule in the coordinate direction in which the 
derivatives are taken, and the trapezoidal rule in the remaining coordinate directions. 

We now have the following convergence result, which follows immediately from the above analysis, to¬ 
gether with the convergence results for the error in approximating p by pf^^ proven in [H [23] . 

Theorem 3.2. Let the permeability k be as in Example\^ or0 and suppose f € Then for any linear 

functional Q bounded on Hq(D), we have 

E[{g{p)-g{pDr]^^‘^ < ch\ 


for any s < 1/2 and q < oo. 

The convergence rate s in Theorem [32] depends on the spatial (Sobolev) regularity of the solution p, and 
is, in the case of discontinuous coefficients k limited above by 1/2. The regularity results proved in [23] only 
considered coefficients that are piecewise continuous with respect to a deterministic partitioning 
of D, but the same arguments can be applied to conclude on the regularity of p in the case of a random 
partitioning. If we assume that for each realisation, the m subdomains are disjoint and convex polygonal, 
and, in the case of zero Dirichlet conditions, no more than 2 subdomains meet at the boundary dD and 
no more than three subdomains meet in the interior, we obtain the optimal convergence rate s < 1/2 in 
Theorem 13.21 

From Theorem l3.21 it now immediately follows from the triangle and reverse triangle inequalities, together 
with V[X] < E[X^] for any random variable X, that assumptions {Ml ) and {M2) are satisfied with a < 1/2 
and /3 < 1, respectively. 

Remark 3.3. / Convergence rates for exact finite element solution) In the presence of quadrature, the rate in 
Theorem l3.2l is optimal. However, the convergence rate of the exact finite element error E[{g{p)—g{p))^))'^]^^‘^ 
is twice that of Theorem 13.21 with the error in functionals g converging with rate 2s, for any s < 1/2. This 
faster convergence rate additionally requires g to be bounded on L'^{D). We will see in section 13131 that we 
in practice often observe this faster convergence rate even with quadrature. 
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Remark 3.4. (Nonlinear functionals) The convergence result in Theorem 13.21 can be shown to hold also for 
non-linear functionals that are Frechet differentiable, under suitable assumptions on the Frechet derivative. 
For more details, see [53]. 

Remark 3.5. (Boundary conditions) The incorporation of non-zero Dirichlet and Neumann conditions is 
handled differently in finite element and finite volume methods. A result of the form in Theorem 13.21 for 
general boundary conditions can again be obtained, by applying the analysis above in interior boxes Bij, 
together with an argument bounding the error over the boundary boxes Bi j (see for example jd] Section 
3.2]). 

3.4 Numerical examples 

In the following two sections, we present numerical results to verify the analysis provided in earlier sections, 
and to show that in fact all three assumptions of Theorem 13.II are satisfied for the model problem considered 
in this paper. In the finite volume discretisation, we use the harmonic average of the permeability values in 
two neighbouring cells to approximate the permeability on the boundary. 

We consider a simple model of a randomly layered medium in two and three spatial dimensions, char¬ 
acterised by three different layers. First consider the two dimensional case, which is illustrated in Figure 
[TJ To sample a realisation of the random medium for a given ui, we draw four uniform random variables 
yi 17(0.8,0.9), j /2 17(0.6, 0.7), ys 17(0.2,0.3) and j /4 17(0.4, 0.5), and then draw straight lines between 

the points (0,j/i) and (l,j/ 2 ), and ( 0 , 2 / 3 ) and ( 1 , 2 / 4 ), respectively. Note that, as shown in Figured] for our 
particular choice of parametrisation, the subdomains created this way are always convex. For 3D simulations, 
we extrude 2D straight lines along the a; 3 -axis to create a 3D layered medium. 



Figure 1: A realisation of the random layers with four uniform random variables 2/1 17(0.8, 0.9), 2/2 ~ 

17(0.6,0.7), 2/3 ~ 17(0.2,0.3) and 2/4 ~ 17(0.4,0.5). 

We consider two different model problems. 

Model Problem 1 (point evaluation of pressure): Firstly, we consider a point value of the pressure. 
Since in higher dimensional space with d = 2,3, point evaluation is not a bounded functional on H^(D), we 
regularise this type of functional by approximating the point value by a local average, 

:= r f udx 

where D* is a small subdomain of D that contains x * (im). As governing equation, we consider the PDF 

—V • (A;(x, w)Vp(x, w)) = 1, for X € (0, l)'^ 
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Figure 2: Results for piecewise constant permeability on random subdomains in 2D. Left-top plot: 
\E[M^^'> (uh*) — M^^\uh)]\ versus 1/h for Model Problem 1. Right-top plot: (uh*) — M^'^'>(uh)]\ 

for Model Problem 2. Left-bottom plot: Y[M^^\uh) — {u 2 h)] versus 1/h for Model Problem 1. Right- 

bottom plot: Y[M^‘^\uh) — M^^'>{u 2 h)] versus 1/h for Model Problem 2. The gradient of the dotted (resp. 
dashed) line is -2 (resp. -1). 


with homogeneous Dirichlet boundary conditions. 

Model Problem 2 (outflow through boundary): Secondly, we consider the PDF 

V • (fc(x, w)Vp(x, w)) = 0, for X e (0, l)'^ (3-11) 

with mixed boundary conditions p\xi=Q = l,p|xi=i = 0 and zero Neumann conditions on the remainder of 
the boundary. We take as quantity of interest the outflow through the boundary xi = 1 

= — / (^(x)fc(x,a;)Vu(x, w) • nds 

JdD 

where </) is a weighting function such that 4>\{xi=i} = 1 and (j)\{xi=o} = 0- Notice that is indeed 

equal to the outflow over the boundary xi = 1 because of the Neumann boundary conditions imposed on p. 

To estimate the errors, we approximate the exact solution it by a reference solution Uh* on a grid with 
mesh size h* = 1/512 in 2D and h* = 1/128 in 3D. 

Figures [2] and [3] show results for the piecewise constant permeability model in Example 1. As described 
previously, we divide the medium into three (random) horizontal layers. For each sample of the layers, we then 
sample from three independent, standard normal random variables Zi,Z 2 and z^, and set the permeability 
values in the three regions to exp(2;i),exp( 2 ; 2 ) and exp( 2 ; 3 ), respectively. 

In Figured! the two plots in the top row show quadratic convergence in h for the error in expected value 
of and linear convergence for M^‘^\ Theorem 13.21 and Remark |3.3I suggest square-root convergence in 
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Figure 3: Results for piecewise constant permeability on random subdomains in 3D. Left-top plot: 
\E[M^^'> (uh*) — M^^\uh)]\ versus 1/h for Model Problem 1. Right-top plot: {uh*) — M^'^'>(uh)]\ 

for Model Problem 2. Left-bottom plot: Y[M^^\uh) — {u 2 h)] versus 1/h for Model Problem 1. Right- 

bottom plot: Y[M^‘^\uh) — M^^'>{u 2 h)] versus 1/h for Model Problem 2. The gradient of the dotted (resp. 
dashed) line is -2 (resp. - 1 ). 


the presence of quadrature, and linear convergence in the absence of quadrature, and so these convergence 
rates are faster than expected. In the two bottom plots, we observe faster than quadratic convergence for 
Y[M^^\uh) — M^^\u 2 h)], and quadratic convergence for Y[M^‘^\uh) — M^‘^\u 2 h)]- In Figure [H we observe 
the same convergence rates in three spatial dimensions. 

Figured shows results for the piecewise correlated field model from Example 2 in two spatial dimensions. 
As before, we divide D = (0,1)^ into three (random) horizontal layers, and model the permeability in the 3 
layers by 2 different log-normal distributions. The parameters in the top and bottom layer are taken to be 
= 0, Ai = 0.3 and af = 1, and for the middle layer we take ^2 = 4, A 2 = 0.1 and = 1) assuming no 
correlation across layers. We use the 2-norm exponential covariance function (12.3p . To generate realisations 
of a log-normal random field, we use the circulant embedding technique, which is an exact and fast method 
of generating samples from stationary Gaussian random fields on a regular grid [261113) . The complexity of 
this sampling method is 0{hj'^ log/i^'^). 

In FigureUl we observe 0{h '^) convergence for the error in expected value for and linear convergence 
for 

3.5 Scalable linear solver (7 = d) in 3D 

In this section, we investigate the scalability of the iterative linear solver, which determines the value of 
7 in Theorem 13.11 Since it is the most challenging, we only consider the case of three spatial dimensions. 
Results in one and two spatial dimensions are similar. As already mentioned in section 13.21 the value of 7 
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Figure 4: Results for piecewise correlated field permeability on random subdomains in 2D. Left plot: 
\E[M^^\uh*) — M^^''{uh)]\ versus 1/h for Model Problem 1. Right plot: \E[M^‘^\uh*) — M^‘^\uh)]\ for 
Model Problem 2. The gradient of the dotted (resp. dashed) line is -3/2 (resp. -1). 


also depends on the sampling method used. The Circulant Embedding method used in the previous section 
has log-linear complexity in the number of degrees of freedom, log Out of the sampling and the 

linear solver, the linear solver usually is the more costly in terms of absolute CPU time. 

The finite volume discretisation of a realisation of model problem (12.21) results in a symmetric positive 
definite (SPD) stiffness matrix. As a linear solver for SPD systems, we use a preconditioned conjugate 
gradient (PCG) method, with the preconditioner given by one cycle of an algebraic multigrid (AMG) solver. 

In 3D simulations, memory requirements for storing the stiffness matrices increases rapidly with the mesh 
width hi, since the number of degrees of freedom in the finite volume method grows like hj'^. Thus, the 
efficiency of the linear solver strongly depends on the number of non-zeros in the matrices. In order to 
achieve optimal computational complexity, we use the aggressive coarsening technique described in |25] . The 
outer PCG iteration is run until the stopping criterion 

||r^||/||ro|| < 10-1° 

is satisfied, where and rg are the residuals obtained at iterations m and 0, respectively. 

Figure [5] shows the average CPU time (in seconds) for the solution of 1000 linear systems using AMG- 
PCG. The systems were generated by discretising the partial differential equation (I3.11|) as in Model Problem 
2, with random coefficients given by a log-normal random field with 2-norm exponential covariance function 
(12.31) with = 1 and A = 0.3. The CPU time increases linearly with the degrees of freedom. Hence, 
AMG-PCG is a scalable linear solver for 3D problems with random coefficients. 

4 Variance reduction through Coarse Grid Variates 

Different variance reduction techniques, such as antithetic variates, control variates and importance sampling, 
have been developed to increase the accuracy of and speed up the simulation process [2] . In this section, we 
introduce a new variance reduction technique specifically designed for multilevel Monte Carlo estimators, 
which we term Coarse Grid Variates (CGV). In particular, the variance reduction technique is designed to 
reduce the variance of the Monte Garlo estimators of the differences Qi — Qi-i, for ^ = 1,... ,L. We will 
describe the method in the context of finite volume discretisations of the Darcy flow equation as discussed 
in earlier sections, but it is easily adaptable to other applications. 
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Figure 5: Average CPU time (in seconds) for the solution of 1000 linear systems with random coefficients, 
solved by using AMG-PCG. 
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Figure 6: ID uniform grid {xq, ... ,Xm^+i} G D = [0,1], with = Xi — Xi-\. Grey circles represent the 
locations of sampling points ke(i). 


4.1 Description of method 

For ease of presentation, let us for the moment assume that the finite volume meshes used in the simulation 
are regular in the sense that nodes are uniformly distributed in each coordinate direction. We also assume 
that the number of cells in each coordinate direction is a power of 2, with s = 2 in the growth condition (13.41) . 
Under these assumptions, the hierarchy of finite volume meshes is nested. We will crucially also assume that 
the distribution of k is stationary in the sense that it is invariant under shifts in d-dimensional space. 

To compute a sample Qi^i — we need to produce a sample of the coefficient k at the cell centres 

of the finite volume meshes on level i and I — Let us denote these coarse and fine samples by ki^i and 
respectively. The aim is now to extract a sample kg-i^i from a given sample k^^i, using the fact that 
the distribution k is stationary. 

We first illustrate our general idea in the one-dimensional setting. Figure [S] shows the computational 
domain D = (0,1), divided into mi cells. The fine grid sample ki^i consists of mi values. The coarse grid 
sample should consist of mi-i = m^/2 values, and since the distribution of k is stationary, it follows 

that the two subvectors of ki^i given by its even and odd entries, respectively, are both acceptable as samples 
ki-i^i. With := [ki{l) ,ki{3 ); • • • , ki{mi - 1)], and kf'^ = [ki{2) ,ki{4) , • • • ,ki{mi)\, both Qf,* 

and Qi^i — Q^ili ^ are valid samples of Qi — Qi-\. Furthermore, we notice that Qi^i — , -I- j)/2 is 

also a valid sample, and this is what we will use as a basis for our variance reduction technique. 

In general d dimensional space, the random field ki on level i has 2^^ sub-vectors k^/\ for j = 1,..., 2^^, 
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Figure 7: 


Variance plot of Yu for Model Problem 2 in 2D,with continuous log-normal permeability. 


which have the distribution required for the random field on level £ — 1. Given a sample of we 
then define the averaged coarse grid quantity 


= ( 4 . 1 ) 

1=1 

and finally the estimator 

E (4-2) 

i—1 

for Yi = Qi — Qi-i- Note that this is still an unbiased estimator of E[lf]. As we will see in sectionthe 
estimator (1121) has a considerably smaller variance than the standard MC estimators for — Qi-\ based 
on only one sample This is due to the fact that the CGV estimator (j4.2p incorporates much more 

information due to the averaging of samples on the coarse level. 

The computational cost of the multilevel estimator based on a telescoping sum of estimators (1121) will 
be less than twice that of of the standard MLMG estimator presented in section 13.21 For each difference 
Qt. — Qi-\i the total cost of the additional computations required is no more than the cost of a single solve 
on the hne level since we have 7 > d. In practice, the sampling cost of the Coarse Grid Variates is less 
than that of 2 Hne-grid samples because no random field needs to be generated on the coarse grid. Instead, 
the random vector already generated on the fine grid is used on the coarse grid. 

4.2 Numerical examples 

In this section we provide numerical results demonstrating the effectiveness of the CGV method. We consider 
the mixed boundary value problem Model Problem 2 already considered in section 13.41 with quantity of 
interest the outflow functional We model fc as a continuous log-normal random field, such that log A: 

has 1-norm exponential covariance function (12.31) with A = 0.3 and = 1. 

Figure 0 shows results for the model problem in two spatial dimensions. In particular, it shows the 
variance of Ye with and without the use of Coarse Grid Variates. The dashed line with diamonds indicates 
the results using CGV. We see that the use of CGV leads to a significantly lower variance, both in terms 
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Figure 8: Variance plot of Yi for Model Problem 2 in 3D,with continuous log-normal permeability. 


of the constant and even a faster decay rate with mesh width hi. Already on the coarsest grid of width 
hi = 1/8, the variance of Yi is reduced by a factor of almost 200. These results suggest that the CGV can 
further lower the rate of growth of the computational cost of the MLMC algorithm. 

In Figure [U we observe a very similar behaviour of Yi for the three dimensional model problem. On the 
coarsest grid hi = 1/8, the variance of Yi is reduced by a factor of 100 when using CGV. We again also 
observe a faster decay rate of Yi with hi when using CGV. 

In numerical simulations with other quantities of interest, such as the point evaluation considered in 
Model Problem 1 from section 1331 the benefits of using the CGV approach were again clear, although not 
as pronounced as in Figures [7] and El The variance reduction attainable from the CGV approach seems to 
depend on the specific model problem and quantity of interest. However, in our tests, the CGV approach 
always resulted in gains over standard MLMC. 


5 Conclusions 

In this work, we considered the application of multilevel Monte Carlo methods to elliptic PDFs with log¬ 
normal random coefficients. We extended the existing theory to cover simple finite volume discretisations of 
the governing equations, and model problems where the location of different layers in the subsurface is also 
subject to uncertainty. Theoretical results were confirmed by numerical simulations. Finally, we proposed 
a new variance reduction technique, termed Coarse Grid Variates, designed for the multilevel Monte Carlo 
method. Numerical simulations show a variance reduction around 2 orders of magnitude compared to 
standard multilevel Monte Carlo for model problems in two and three spatial dimensions. 
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